#!/usr/bin/python
# -*- coding: utf-8 -*-
__author__="Boris Avdeev"
__date__ ="$Mar 8, 2011 10:09:48 AM$"

#TODO: Migrate into AFTdata class?



import sys
import numpy as np
import matplotlib.mlab as mlab
import matplotlib.pyplot as plt

def est_dpar(dpars,sample):
    """Returns mean and standard deviation and plots histogram of values with a normal pdf on top for diagnostic"""
    Dpar=np.mean(dpars)
    Dpar_std=np.std(dpars)
    plt.hist(dpars, normed=True, facecolor='gray', alpha=0.75)
    x = np.linspace(min(dpars),max(dpars))
    y = mlab.normpdf(x, Dpar, Dpar_std)
    plt.plot(x, y, 'b-', linewidth=1)
    plt.xlabel('Dpar')
    plt.ylabel('Probability')
    plt.title('Sample %s, N = %i' % (sample,len(dpars)))
    #plt.axis([40, 160, 0, 0.03])
    plt.grid(True)
    plt.savefig("%s.pdf"%sample)
    return Dpar, Dpar_std

if __name__ == '__main__':
  """ Create a qtqt config file from AFTSolve .A2G and . AGE files and AHe .csv spreadsheet.
      Code assumes many things about input format and aslo hard-codes some parameters.
      Apart from format conversion, it computes induced track densities (required by QTQt)
      from LA-ICP-MS P value (P = U/Ca, Donelick et al., 2005) """
  aft = False
  ns, lens = [], []
  zeta, rho_d, n_d, Dpar, Dpar_std = 0, 0, 0, 0, 0
  sample = sys.argv[1]
  try:
    cnt_file = open("A%sA.AGE" % sample)
    len_file = open("L%sA.A2G" % sample)
    aft = True
  except:
    print "No AFT data for sample %s" % sample
  if aft:
    cnt_head = cnt_file.readlines(6)
    len_head = len_file.readlines(3)
    cnt_file.seek(0)
    len_file.seek(0)
    sample = cnt_head[1].strip()
    if(sample!=len_head[0].strip()):
      message = "Mismatch in sample names between %s and %s files" % (cnt_file.name, len_file.name)
      raise Exception, message
    zeta_ms = float(cnt_head[4].split()[0])  # zeta_ms
    zeta  = zeta_ms * 20  # Number to approximately match real value.  Should not matter without error handling.
    rho_d = 1000000       # Number to approximately match real counts. Should not matter without error handling.
    n_d = 1000            # Number to approximately match real counts. Should not matter without error handling.
    counts = np.genfromtxt(cnt_file,skiprows=6)

    Dpar, Dpar_std = est_dpar(counts[:,4], sample)

    ns = counts[:,0]
    #          area  *  P  * rho_d * zeta / zeta_ms
    ni = counts[:,1] * counts[:,2] * rho_d * zeta / zeta_ms 

    lengths= np.genfromtxt(len_file,skiprows=3)
    lens = lengths[:,0]
    angs = lengths[:,1]

  n_tT = 0  # Put constraints in the total profile
  auth = np.genfromtxt("auth_sg.csv", dtype=None, delimiter=',',skiprows=0,names=True)
  auth = auth[ auth['sample']==sample.lower() ]
  he_ncc = 22413.6 * auth['He']

  
  out_file = open("%s.qtqt" % sample,"w")
  out_file.write('%s\n'%sample)
  out_file.write("0 0 0\n") # x y z
  out_file.write("%i %i %i %f %f %i\n" % (n_tT, len(lens), len(ns), zeta, rho_d, n_d))
  out_file.write("105\n") # Ketcham2007 model
  out_file.write("0 %f %f\n" % (Dpar, Dpar_std))
  out_file.write("1 16.3\n") 
  out_file.write("1\n")  # Project lengths
  out_file.write("1\n")  # Cf
  out_file.write("0\n")  # Donelick's etching
  #out_file.write("10 -10\n") # Surface T
  out_file.write("0 0\n") # FT age
  out_file.write("0 0\n") # TL mean
  out_file.write("0 0\n") # TL std.dev.
  for i in range(len(ns)):
    out_file.write("%f %f\n" % (ns[i], ni[i])) 
    #out_file.write("%i %i\n" % (ns[i], ni[i]+0.5))  # replace with %f %f and remove 0.5 when Kerry fixes his code
  for i in range(len(lens)):
    out_file.write("%f %f\n" % (lens[i], angs[i]))
  out_file.write("%i\n" % len(auth)) # AHe age No
  out_file.write("2\n") # Use Flowers rad. damage model
  for i in range(len(auth)):
    out_file.write("%f %f %f %f %f %f %f 0\n" % 
        (he_ncc[i], auth[i]['U'], auth[i]['Th'], auth[i]['raw'], auth[i]['raw']*0.15, auth[i]['l'], auth[i]['r']  ))
    out_file.write("20 0.005 138000\n")
     
  
  
  out_file.close()
  

